perm filename FFT.SAV[3,ALS] blob
sn#201661 filedate 1971-06-20 generic text, type T, neo UTF8
00100 BEGIN "FAST"
00200 COMMENT This program does a series of fast Fourier transforms ON
00250 Segments of data taken from an
00300 input file specified from the terminal
00400 and creates an output file containing a series of
00500 power spectra. The size of the segments to be used
00600 used is specified by typing in a value for M,
00700 which must be ≤ 8;
00800
01000 REQUIRE "DPYSUB.HDR[1,PDQ]" SOURCE_FILE;
01100 DEFINE DPYSIZ="1000";
01200 INTEGER ARRAY DPYBUF[1:DPYSIZ];
01300
01350 EXTERNAL STRING PROCEDURE DATIME;
01375 EXTERNAL PROCEDURE TIMSET;
01400 EXTERNAL STRING PROCEDURE STRIN(STRING S);
01500 EXTERNAL INTEGER PROCEDURE ININT(STRING S);
01525 EXTERNAL REAL PROCEDURE RUNTIM;
01550 EXTERNAL STRING PROCEDURE INCHWL;
01575 EXTERNAL PROCEDURE WAIT(INTEGER I);
01600
01700 FORTRAN REAL PROCEDURE ALOG10(REAL X);
01800 FORTRAN REAL PROCEDURE COS(REAL X);
01900 FORTRAN REAL PROCEDURE SIN(REAL X);
02000 FORTRAN REAL PROCEDURE SQRT(REAL X);
02100
02200 REAL ARRAY A,B,C[0:256];
02250 REAL ARRAY W[0:128];
02300
02400 STRING FILEI,TFILEI,CHAR1,CHAR2,CHAR3,CHAR4,CHAR5,CHAR6,CHAR7,LOGFRE,CEPS;
02500
02600 INTEGER ARRAY LFILE[0:'177];
02700 INTEGER ARRAY D[0:256];
02800 INTEGER ARRAY DATBUF[0:1023];
02900
03000 INTEGER AIFORM,AOFORM,EOF,IEOF,I,J,K,M,TM,N,BPT,SEGC,SEGCNT,SEGTOT,SAMINT;
03050 INTEGER XPOS,YPOS,MAX;
03075 EXTERNAL REAL RTIM,ATIM;
03100
03200 DEFINE DATSIZ="1024";
03300 DEFINE BPS="12";
03350 DEFINE LEFT="24";
03375 DEFINE DI="2↑24";
03400
03500 LABEL START;
03600
00100 PROCEDURE UNPACK(REAL ARRAY A,B;INTEGER ARRAY DATBUF;REFERENCE INTEGER BPT;INTEGER N,SAMINT);
00200 BEGIN
00300 COMMENT This procedure accepts a segment of digitized speech of
00400 length 2↑(m+1). It expands this segment and
00500 loads it into A andB in real form;
00600 INTEGER I,J;
00800 FOR I←0 STEP 1 UNTIL N-1 DO
00815 BEGIN
00826 A[I]←((ILDB(BPT) LSH LEFT)%DI);
00830 IF (SAMINT)=2 THEN J←ILDB(BPT);
00860 END;
00900 FOR I←0 STEP 1 UNTIL N-1 DO
00930 BEGIN
00940 B[I]←((ILDB(BPT) LSH LEFT)%DI);
00950 IF (SAMINT)=2 THEN J←ILDB(BPT);
00960 END;
01060 END "UNPACK";
01100
01200 PROCEDURE DPYDMP;
01300 BEGIN STRING FILE;INTEGER DSIZ,CCCHN;
01400 IF ¬(FILE←STRIN("CAL-COMP FILE=")) THEN RETURN;
01500 OPEN(CCCHN←GETCHAN,"DSK",'14,0,1,0,0,0);
01600 ENTER(CCCHN,FILE,0);
01700 DPYPARS;DSIZ←DPYBUF[2]+3;
01800 ARRYOUT(CCCHN,DPYBUF[1],DSIZ);
01900 RELEASE(CCCHN);
02000 END "DPYDMP";
02100
02200
02300 PROCEDURE FFT(REAL ARRAY A,B;INTEGER N,M,KS);
02400 BEGIN
02500 COMMENT COMPUTES THE FFT FOR ONE VARIABLE OF DIMENSION 2↑M;
02600 INTEGER K0,K1,K2,K3,SPAN,J,JJ,K,KB,KN,MM,MK;
02700 REAL RAD,C1,C2,C3,S1,S2,S3,CK,SK,SQ;
02800 REAL A0,A1,A2,A3,B0,B1,B2,B3;
02900 INTEGER ARRAY C[0:M];
03000 LABEL L,L2,L3,L4,L5,L6;
03100 SQ←0.707106781187;
03200 SK←0.382683432366;
03300 CK←0.92387953251;
03400 C[M]←KS; MM←(M%2)*2; KN←0;
03500 FOR K←M-1 STEP -1 UNTIL 0 DO C[K]←C[K+1]/2;
03600 RAD←6.28318530718/(C[0]*KS); MK←M-5;
03700 L: KB←KN; KN←KN+KS;
03800 IF MM≠N THEN
03900 BEGIN
04000 K2←KN; K0←C[MM]+KB;
04100 L2: K2←K2-1; K0←K0-1;
04200 A0←A[K2]; B0←B[K2];
04300 A[K2]←A[K0]-A0; A[K0]←A[K0]+A0;
04400 B[K2]←B[K0]-B0; B[K0]←B[K0]+B0;
04500 IF K0>KB THEN GO TO L2;
04600 END;
04700 C1←1.0; S1←0;
04800 JJ←0; K←MM-2; J←3;
04900 IF K≥0 THEN GO TO L4 ELSE GO TO L6;
05000 L3: IF C[J]≤JJ THEN
05100 BEGIN
05200 JJ←JJ-C[J]; J←J-1;
05300 IF C[J]≤JJ THEN
05400 BEGIN
05500 JJ←JJ-C[J]; J←J-1; K←K+2;
05600 GO TO L3;
05700 END
05800 END;
05900 JJ←C[J]+JJ; J←3;
06000 L4: SPAN←C[K];
06100 IF JJ≠0 THEN
06200 BEGIN
06300 C2←JJ*SPAN*RAD; C1←COS(C2); S1←SIN(C2);
06400 L5: C2←C1↑2-S1↑2; S2←2.0*C1*S1;
06500 C3←C2*C1-S2*S1; S3←C2*S1+S2*C1;
06600 END;
06700 FOR K0←KB+SPAN-1 STEP -1 UNTIL KB DO
06800 BEGIN
06900 K1←K0+SPAN; K2←K1+SPAN; K3←K2+SPAN;
07000 A0←A[K0]; B0←B[K0];
07100 IF S1=0 THEN
07200 BEGIN
07300 A1←A[K1]; B1←B[K1];
07400 A2←A[K2]; B2←B[K2];
07500 A3←A[K3]; B3←B[K3];
07600 END
07700 ELSE
07800 BEGIN
07900 A1←A[K1]*C1-B[K1]*S1;
08000 B1←A[K1]*S1+B[K1]*C1;
08100 A2←A[K2]*C2-B[K2]*S2;
08200 B2←A[K2]*S2+B[K2]*C2;
08300 A3←A[K3]*C3-B[K3]*S3;
08400 B3←A[K3]*S3+B[K3]*C3;
08500 END;
08600 A[K0]←A0+A2+A1+A3; B[K0]←B0+B2+B1+B3;
08700 A[K1]←A0+A2-A1-A3; B[K1]←B0+B2-B1-B3;
08800 A[K2]←A0-A2-B1+B3; B[K2]←B0-B2+A1-A3;
08900 A[K3]←A0-A2+B1-B3; B[K3]←B0-B2-A1+A3;
09000 END;
09100 IF K>0 THEN BEGIN K←K-2; GO TO L4; END;
09200 KB←K3+SPAN;
09300 IF KB<KN THEN
09400 BEGIN
09500 IF J=0 THEN BEGIN K←2; J←MK; GO TO L3; END;
09600 J←J-1; C2←C1;
09700 IF J=1 THEN
09800 BEGIN C1←C1*CK+S1*SK; S1←S1*CK-C2*SK; END
09900 ELSE BEGIN C1←(C1-S1)*SQ; S1←(C2+S1)*SQ; END;
10000 GO TO L5;
10100 END;
10200 L6: IF KN<N THEN GO TO L;
10300 END "FFT";
10400
10500
10600
10700
10800 PROCEDURE REVFFT(REAL ARRAY A,B;INTEGER N,M,KS);
10900 BEGIN
11000 COMMENT COMPUTES THE FFT FOR ONE VARIABLE OF DIMENSION 2↑M IN A
11100 MULTIVARIATE TRANSFORM.
11200 IF N=2↑M AND K=1 THEN A SINGLE-VARIATE TRANSFORM IS COMPUTED;
11300 INTEGER K0,K1,K2,K3,K4,SPAN,NN,J,JJ,K,KB,NT,KN,MK;
11400 REAL RAD,C1,C2,C3,S1,S2,S3,CK,SK,SQ;
11500 REAL A0,A1,A2,A3,B0,B1,B2,B3,RE,IM;
11600 INTEGER ARRAY C[0:M];
11700 LABEL L,L2,L3,L4,L5,L6;
11800 SQ←0.707106781187;
11900 SK←0.382683432366;
12000 CK←0.92387953251;
12100 C[0]←KS; KN←0; K4←4*KS; MK←M-4;
12200 FOR K←1 STEP 1 UNTIL M DO C[K]←KS←KS+KS;
12300 RAD←3.1415926536/(C[0]*KS);
12400 L: KB←KN+K4; KN←KN+KS;
12500 IF M=1 THEN GO TO L5;
12600 K←JJ←0; J←MK; NT←3;
12700 C1←1.0; S1←0;
12800 L2: SPAN←C[K];
12900 IF JJ≠0 THEN
13000 BEGIN
13100 C2←JJ*SPAN*RAD; C1←COS(C2); S1←SIN(C2);
13200 L3: C2←C1↑2-S1↑2; S2←2*C1*S1;
13300 C3←C2*C1-S2*S1; S3←C2*S1+S2*C1;
13400 END
13500 ELSE S1←0;
13600 K3←KB-SPAN;
13700 L4: K2←K3-SPAN; K1←K2-SPAN; K0←K1-SPAN;
13800 A0←A[K0]; B0←B[K0];
13900 A1←A[K1]; B1←B[K1];
14000 A2←A[K2]; B2←B[K2];
14100 A3←A[K3]; B3←B[K3];
14200 A[K0]←A0+A1+A2+A3; B[K0]←B0+B1+B2+B3;
14300 IF S1=0 THEN
14400 BEGIN
14500 A[K1]←A0-A1-B2+B3; B[K1]←B0-B1+A2-A3;
14600 A[K2]←A0+A1-A2-A3; B[K2]←B0+B1-B2-B3;
14700 A[K3]←A0-A1+B2-B3; B[K3]←B0-B1-A2+A3;
14800 END
14900 ELSE
15000 BEGIN
15100 RE←A0-A1-B2+B3; IM←B0-B1+A2-A3;
15200 A[K1]←RE*C1-IM*S1; B[K1]←RE*S1+IM*C1;
15300 RE←A0+A1-A2-A3; IM←B0+B1-B2-B3;
15400 A[K2]←RE*C2-IM*S2; B[K2]←RE*S2+IM*C2;
15500 RE←A0-A1+B2-B3; IM←B0-B1-A2+A3;
15600 A[K3]←RE*C3-IM*S3; B[K3]←RE*S3+IM*C3;
15700 END;
15800 K3←K3+1; IF K3<KB THEN GO TO L4;
15900 NT←NT-1;
16000 IF NT≥0 THEN
16100 BEGIN
16200 C2←C1;
16300 IF NT=1 THEN
16400 BEGIN C1←C1*CK+S1*SK; S1←S1*CK-C2*SK; END
16500 ELSE BEGIN C1←(C1-S1)*SQ; S1←(C2+S1)*SQ; END;
16600 KB←KB+K4; IF KB≤KN THEN GO TO L3 ELSE GO TO L5;
16700 END;
16800 IF NT=-1 THEN BEGIN K←2; GO TO L2; END;
16900 IF C[J]≤JJ THEN
17000 BEGIN
17100 JJ←JJ-C[J]; J←J-1;
17200 IF C[J]≤JJ THEN
17300 BEGIN JJ←JJ-C[J]; J←J-1; K←K+2; END
17400 ELSE BEGIN JJ←C[J]+JJ; J←MK;END;
17500 END
17600 ELSE BEGIN JJ←C[J]+JJ; J←MK; END;
17700 IF J<MK THEN GO TO L2; K←0; NT←3;
17800 KB←KB+K4; IF KB≤KN THEN GO TO L2;
17900 L5: K←(M%2)*2;
18000 IF K≠M THEN
18100 BEGIN
18200 K2←KN; K0←J←KN-C[K];
18300 L6: K2←K2-1; K0←K0-1;
18400 A0←A[K2]; B0←B[K2];
18500 A[K2]←A[K0]-A0; A[K0]←A[K0]+A0;
18600 B[K2]←B[K0]-B0; B[K0]←B[K0]+B0;
18700 IF K2>J THEN GO TO L6;
18800 END;
18900 IF KN<N THEN GO TO L;
19000 END "REVFFT";
19100
19200 PROCEDURE REORDER(REAL ARRAY A,B;INTEGER N,M,KS,REEL);
19300 BEGIN
19400 COMMENT PERMUTES DATA FROM NORMAL TO REVERSE BINARY ORDER AND BACK;
19500 INTEGER I,J,JJ,K,KK,KB,K2,KU,LIM,P;
19600 REAL T;
19700 INTEGER ARRAY C,LST[0:M];
19800 LABEL L,L2,L3,L4;
19900 C[M]←KS;
20000 FOR K←M STEP -1 UNTIL 1 DO C[K-1]←C[K]%2;
20100 P←J←M-1; I←KB←0;
20200 IF REEL THEN
20300 BEGIN
20400 KU←N-2;
20500 FOR K←0 STEP 2 UNTIL KU DO
20600 BEGIN T←A[K+1]; A[K+1]←B[K]; B[K]←T; END;
20700 END
20800 ELSE M←M-1;
20900 LIM←(M+2)%2; IF P≤0 THEN GO TO L4;
21000 L: KU←K2←C[J]+KB; JJ←C[M-J]; KK←KB+JJ;
21100 L2: K←KK+JJ;
21200 L3: T←A[KK]; A[KK]←A[K2]; A[K2]←T;
21300 T←B[KK]; B[KK]←B[K2]; B[K2]←T;
21400 KK←KK+1; K2←K2+1;
21500 IF KK<K THEN GO TO L3;
21600 KK←KK+JJ; K2←K2+JJ;
21700 IF KK<KU THEN GO TO L2;
21800 IF J>LIM THEN
21900 BEGIN
22000 J←J-1; I←I+1;
22100 LST[I]←J; GO TO L;
22200 END;
22300 KB←K2;
22400 IF I>0 THEN
22500 BEGIN J←LST[I]; I←I-1; GO TO L; END;
22600 IF KB<N THEN BEGIN J←P; GO TO L; END;
22700 L4: ;
22800 END "REORDER";
22900
23000
23100 PROCEDURE RTRAN(REAL ARRAY A,B;INTEGER N,EVALUATE);
23200 BEGIN
23300 COMMENT IF EVALUATE IS FALSE THIS PROCEDURE UNSCRAMBLES THE SINGLE VARIATE
23400 COMPLEX TRANSFORM ;
23500 INTEGER K,NK,NH;
23600 REAL AA,AB,BA,BB,RE,IM,CK,SK,DC,DS,R;
23700 NH←N%2; R←3.1415926536/N;
23800 DS←SIN(R); R←-(2*SIN(0.5*R))↑2;
23900 DC←-0.5*R; CK←1.0; SK←0;
24000 IF EVALUATE THEN
24100 BEGIN
24200 CK←-1.0; DC←-DC;
24300 END
24400 ELSE
24500 BEGIN
24600 A[N]←A[0]; B[N]←B[0];
24700 END;
24800 FOR K←0 STEP 1 UNTIL NH DO
24900 BEGIN
25000 NK←N-K;
25100 AA←A[K]+A[NK]; AB←A[K]-A[NK];
25200 BA←B[K]+B[NK]; BB←B[K]-B[NK];
25300 RE←CK*BA+SK*AB; IM←SK*BA-CK*AB;
25400 B[NK]←IM-BB; B[K]←IM+BB;
25500 A[NK]←AA-RE; A[K]←AA+RE;
25600 DC←R*CK+DC; CK←CK+DC;
25700 DS←R*SK+DS; SK←SK+DS;
25800 END;
25900 END "RTRAN";
26000
26100
26200 PROCEDURE RFOUR(REAL ARRAY A,B;INTEGER M,INVERSE);
26300 BEGIN
26400 COMMENT COMPUTES THE FFT OF 2↑(M+1) REAL DATA POINTS;
26500 INTEGER N,J;
26600 REAL P;
26700 N←2↑M;
26800 IF INVERSE THEN
26900 BEGIN
27000 RTRAN(A,B,N,TRUE);
27100 FOR J←N-1 STEP -1 UNTIL 0 DO
27200 B[J]←-B[J];
27300 FFT(A,B,N,M,N);
27400 FOR J←N-1 STEP -1 UNTIL 0 DO
27500 BEGIN A[J]←0.5*A[J]; B[J]←-0.5*B[J] END;
27600 REORDER(A,B,N,M,N,TRUE);
27700 END
27800 ELSE
27900 BEGIN
28000 REORDER(A,B,N,M,N,TRUE);
28100 REVFFT(A,B,N,M,1); P←0.5/N;
28200 FOR J←N-1 STEP -1 UNTIL 0 DO
28300 BEGIN A[J]←P*A[J]; B[J]←P*B[J] END;
28400 RTRAN(A,B,N,FALSE);
28500 END;
28600 END "RFOUR";
28700
28800
28900 PROCEDURE CFOUR(REAL ARRAY A,B;INTEGER M,INVERSE);
29000 BEGIN
29100 COMMENT COMPUTES THE FFT OF 2↑M COMPLEX DATA VALUES Xi IF INVERSE IS TRUE;
29200 INTEGER N,J;
29300 REAL P,Q;
29400 N←2↑M; P←Q←1.0/SQRT(N);
29500 IF INVERSE THEN
29600 BEGIN
29700 Q←-Q;
29800 FOR J←N-1 STEP -1 UNTIL 0 DO B[J]←-B[J];
29900 END;
30000 FFT(A,B,N,M,N); REORDER(A,B,N,M,N,FALSE);
30100 FOR J←N-1 STEP -1 UNTIL 0 DO
30200 BEGIN A[J]←A[J]*P; B[J]←B[J]*Q; END;
30300 END "CFOUR";
30400
30500 PROCEDURE POWER(REAL ARRAY A,B,C;INTEGER N);
30600 BEGIN
30700 COMMENT THIS COMPUTES THE POWER SPECTRUM OF THE SIN AND COS SERIES IN A,B;
30800 INTEGER I;
30900 FOR I←0 STEP 1 UNTIL N DO
31000 C[I]←SQRT(A[I]↑2 + B[I]↑2);
31100 END "POWER";
31200
31300 PROCEDURE ARRDIS(REAL ARRAY A; INTEGER N,XPOS,YPOS;STRING ID);
31400 BEGIN
31500 COMMENT DISPLAYS A HISTOGRAM OF THE FIRST N VALUES OF ARRAY A AT XPOS,YPOS;
31600 INTEGER I,J,SP;
31700 INTEGER LY,DY;
31800 REAL MAX;
31900 MAX←0;
32000 FOR I←0 STEP 1 UNTIL N DO
32100 IF ABS(A[I])>MAX THEN MAX←ABS(A[I]);
32200 MAX←MAX/360;
32300 SP←1024%N; COMMENT HORIZONTAL SPACING;
32400 AIVECT(XPOS,YPOS); RVECT(1023,0); RIVECT(-1023,0); RVECT(0,360); RIVECT(0,-360);
32500 LY←A[0]/MAX+YPOS;
32600 AIVECT(XPOS,LY);
32700 FOR I←1 STEP 1 UNTIL N DO
32800 BEGIN
32900 DY←A[I]/MAX+YPOS-LY;
33000 LY←LY+DY;
33100 RVECT(SP,DY);
33200 END;
33300 AIVECT(XPOS,YPOS);
33400 FOR I←1 STEP 1 UNTIL 10 DO
33500 BEGIN
33520 RVECT(0,-15); COMMENT INSERT HORIZONTAL SCALE;
33540 RIVECT(26,15);
33560 RVECT(0,-5);
33580 RIVECT(26,5);
33600 RVECT(0,-10);
33700 RIVECT(26,10);
33750 RVECT(0,-5);
33775 RIVECT(26,5);
33800 END;
33850 RVECT(0,-15);
33900 AIVECT(XPOS,YPOS-40);
33940 DPYSST("0 1 2 3 4 5 6 7 8 9 10");
33980 AIVECT(XPOS,YPOS-60);
34000 DPYSST(ID);
34100 END "ARRDIS";
34150
34200 PROCEDURE DUBDIS(REAL ARRAY A,B; INTEGER N,XPOS,YPOS;STRING ID);
34300 BEGIN
34400 COMMENT DISPLAYS A HISTOGRAM OF THE FIRST N VALUES OF ARRAY A AT XPOS,YPOS;
34500 INTEGER I,J,SP;
34600 INTEGER LY,DY;
35300 MAX←MAX/250;
35400 SP←512%N; COMMENT HORIZONTAL SPACING;
35500 AIVECT(XPOS,YPOS); RVECT(1023,0); RIVECT(-1023,0); RVECT(0,250); RIVECT(0,-250);
35600 LY←A[0]%10+YPOS;
35700 AIVECT(XPOS,LY);
35800 FOR I←1 STEP 1 UNTIL N DO
35900 BEGIN
36000 DY←A[I]%10+YPOS-LY;
36100 LY←LY+DY;
36200 RVECT(SP,DY);
36300 END;
36400 FOR I←1 STEP 1 UNTIL N DO
36500 BEGIN
36600 DY←B[I]/10+YPOS-LY;
36700 LY←LY+DY;
36800 RVECT(SP,DY);
36900 END;
37000 AIVECT(XPOS,YPOS-40);
37100 DPYSST(ID);
37200 END "DUBDIS";
37300
38000 REAL PROCEDURE TIMDPY(REAL XPOS,YPOS,RTIM;STRING TFILE);
38050 BEGIN
38150 REAL X;
38275 X←RUNTIM-RTIM;
38300 IF CHAR1="Y" THEN
38400 BEGIN
38500 AIVECT(XPOS,YPOS);
38600 DPYSST(TFILE&CVS(X));
38700 END
38800 ELSE
38900 OUTSTR(TFILE&CVS(X));
39100 END "TIMDPY";
39200
40000 INTEGER PROCEDURE FLOW(INTEGER ARRAY DATBUF;INTEGER XPOS,YPOS,M,MAX);
40050 BEGIN
40100 INTEGER ARRAY A[0:512];
40150 INTEGER ARRAY B[0:256];
40200 INTEGER I,J,K,L,BPT,EOF,SEGC,AIFORM,SP,LY,DY;
40250 AIFORM←1;
40300 SP←2; COMMENT HORIZONTAL SPACING;
40400 OPEN(AIFORM,"DSK",'10,10,0,0,0,EOF);
40500 LOOKUP(AIFORM,FILEI&".DMD",0);
40600 ARRYIN(AIFORM,DATBUF[0],'200); COMMENT INPUT THE HEADER;
40700 EOF←0; SEGCNT←0; SEGC←0;
40800 MAX←2048%MAX;
40900 WHILE EOF=0 DO
41000 BEGIN
41100 ARRYIN(AIFORM,DATBUF[0],DATSIZ);
41200 IF EOF≠0 THEN
41300 BEGIN
41400 J←EOF LAND '777777;
41500 FOR I←J STEP 1 UNTIL 1023 DO DATBUF[I]←0;
41600 END;
41700 BPT←POINT(BPS,DATBUF[0],-1);
41800 FOR I←0 STEP 1 UNTIL 511 DO A[I]←((ILDB(BPT) LSH LEFT)%DI);
41900 FOR I←1 STEP 1 UNTIL 3*DATSIZ%256-2 DO
42000 BEGIN
42100 FOR J←0 STEP 1 UNTIL 255 DO B[J]←((ILDB(BPT) LSH LEFT)%DI);
42200 FOR K←0 STEP 2↑M UNTIL 255 DO
42300 BEGIN
42400 DPYSET(DPYBUF);
42500 AIVECT(XPOS+512-K*SP,YPOS); RVECT(511,0); RIVECT(-400,-200);
42550 DPYSST("Segment number "&CVS(J));
42600 LY←A[K]/MAX+YPOS;
42700 AIVECT(XPOS,LY);
42800 FOR L←K+1 STEP 1 UNTIL 511 DO
42900 BEGIN
43000 DY←A[L]/MAX+YPOS-LY;
43100 LY←LY+DY;
43200 RVECT(SP,DY);
43300 END;
43400 FOR L←0 STEP 1 UNTIL K DO
43500 BEGIN
43600 DY←B[L]/MAX+YPOS-LY;
43700 LY←LY+DY;
43800 RVECT(SP,DY);
43900 END;
44000 DPYOUT(1);
44100 END;
44200 FOR J←0 STEP 1 UNTIL 255 DO
44300 BEGIN
44400 A[J]←A[J+256];
44500 A[J+256]←B[J];
44600 END;
44700 END;
44800 END;
44900 END"FLOW";
00100 PROCEDURE LOGDIS(REAL ARRAY A; INTEGER ARRAY D; INTEGER N,XPOS,YPOS;STRING ID);
00200 BEGIN
00300 COMMENT DISPLAYS A LOG FREQUENCY HISTOGRAM OF THE FIRST N VALUES OF ARRAY A AT XPOS,YPOS;
00400 INTEGER I,J,SP;
00500 INTEGER LY,DY;
01200 AIVECT(XPOS,YPOS); RVECT(1023,0); RIVECT(-1023,0); RVECT(0,360); RIVECT(0,-360);
01300 LY←YPOS;
01400 AIVECT(XPOS,LY);
01500 FOR I←2 STEP 1 UNTIL N*210%256 DO
01600 BEGIN
01650 SP←D[I]-D[I-1]; COMMENT HORIZONTAL SPACING;
01675 IF A[I]≤1 THEN DY←YPOS-LY ELSE
01700 DY←100*ALOG10(A[I])+YPOS-LY;
01800 LY←LY+DY;
01850 IF I=2 THEN RIVECT(SP,DY) ELSE
01900 RVECT(SP,DY);
02000 END;
02100 AIVECT(XPOS,YPOS);
02150 RVECT(0,-15);
02175 RIVECT(92,15);
02200 FOR I←1 STEP 1 UNTIL 7 DO
02300 BEGIN
02400 RVECT(0,-15); COMMENT INSERT HORIZONTAL SCALE;
02500 RIVECT(33,15);
02600 RVECT(0,-5);
02700 RIVECT(33,5);
02800 RVECT(0,-10);
02900 RIVECT(33,10);
03000 RVECT(0,-5);
03100 RIVECT(34,5);
03200 END;
03300 RVECT(0,-15);
03400 AIVECT(XPOS,YPOS-40);
03450 IF SAMINT=2 THEN
03475 DPYSST("0 32 64 128 256 512 1024 2048 4096")
03487 ELSE
03500 DPYSST("0 64 128 256 512 1024 2048 4096 8192");
03600 AIVECT(XPOS,YPOS-60);
03700 DPYSST(ID);
03800 END "LOGDIS";
03900
04000 PROCEDURE LOGMAG(REAL ARRAY A,B;INTEGER M);
04100 BEGIN
04200 FOR I←0 STEP 1 UNTIL 2↑M DO
04300 BEGIN
04400 A[I]←30*ALOG10(A[I]);
04500 B[I]←30*ALOG10(A[I]);
04600 END;
04700 END "LOGMAG";
04800
04900 PROCEDURE CEPSTRUM(REAL ARRAY A,B;INTEGER M);
05000 BEGIN
05100 RFOUR(A,B,M,0);
05200 LOGMAG(A,B,M);
05300 RFOUR(A,B,M,1);
05400 END "CEPSTRUM";
05500
05510 PROCEDURE WINTAB(REAL ARRAY W;INTEGER N);
05512 BEGIN
05513 INTEGER I;
05514 REAL PI;
05516 PI←3.1315926536;
05518 FOR I←0 STEP 1 UNTIL N%2 DO
05520 W[I]←(1-COS(2*PI*I/N))/2;
05522 END "WINTAB";
05622
06000 PROCEDURE WINDOW(REAL ARRAY A,B;INTEGER N,K);
06100 BEGIN
06200 COMMENT This procedure applies a Tukey Interim Window to arrays A and B
06250 windowing the k data items from each end ;
06300 INTEGER I;
06600 FOR I←0 STEP 1 UNTIL K DO A[I]←A[I]*W[(N-1)*I/(2*K)];
06700 FOR I←N-1-K STEP 1 UNTIL N-1 DO B[I]←B[I]*W[(N-1)*(N-1-I)/(2*K)]
06800 END "WINDOW";
06900
00100 AIFORM←1; AOFORM←2;
00200
00300 START:
00400 IF (TFILEI←STRIN("DATA FILE="))≠NULL THEN FILEI←TFILEI;
00500 CLOSE(AIFORM);
00530 SAMINT←1;
00600 IF (TM←ININT("M="))≠0 THEN M←TM;
00700 IF M>9 THEN GO TO START;
00710 IF M<6 THEN
00720 BEGIN
00730 FLOW(DATBUF,-511,0,M,480);
00740 GO TO START;
00745 END;
00750 IF M=9 THEN
00770 BEGIN
00790 M←M-1; SAMINT←2;
00792 OUTSTR("Taken as M=8 and sampling interval of 2"&'15&'12);
00794 END;
00800
00900 N←2↑M; LFILE[1]←DATSIZ%(6*N);
01000 OPEN(AIFORM,"DSK",'10,10,0,0,0,EOF);
01100 LOOKUP(AIFORM,FILEI&".DMD",0);
01200 ARRYIN(AIFORM,LFILE[0],'200); COMMENT INPUT THE HEADER;
01400 CHAR1←"Y"; CHAR2←"Y"; CHAR3←"N"; CHAR4←"N"; CHAR5←"N";
01500 CHAR6←"Y"; LOGFRE←"Y"; CEPS←"Y";
01600 IF (CHAR7←STRIN("Standard options YorN="))="N" THEN
01700 BEGIN
01800 CHAR1←STRIN("Displays YorN =");
01900 CHAR2←STRIN("FFT wanted YorN =");
01925 IF CHAR2="Y" THEN
01950 LOGFRE←STRIN("Octave scale YorN =");
01975 CEPS←STRIN("CEPSTRUM YorN =");
02000 CHAR3←STRIN("Calcom stop YorN=");
02100 CHAR4←STRIN("Stop anyway YorN=");
02200 CHAR5←STRIN("Save Spectra YorN =");
02300 CHAR6←STRIN("Want timings YorN =");
02400 END;
02410
02420 IF LOGFRE="Y" THEN
02430 BEGIN
02460 FOR I←1 STEP 1 UNTIL 255 DO
02482 D[I]←1024*ALOG10(I*256/N)/ALOG10(210);
02500 D[0]←0; D[1]←0;
02505 END;
02552
02576 WINTAB(W,N); COMMENT PREPARE TABLE FOR WINDQW PROCEDURE;
02588
02600 IF CHAR5="Y" THEN
02700 BEGIN
02800 OPEN(AOFORM,"DSK",'10,0,10,0,0,0);
02900 ENTER (AOFORM,TFILEI&".POW",0);
03000 END;
03100
03200 EOF←0; SEGC←0; SEGCNT←0;
03250 SEGTOT←(LFILE[0]*3+2*N-1)%(2*N);
03300 WHILE EOF=0 DO
03400 BEGIN
03500 ARRYIN(AIFORM,DATBUF[0],DATSIZ);
03600 IF EOF≠0 THEN
03700 BEGIN
03800 J←EOF LAND '777777;
03900 FOR I←J STEP 1 UNTIL N-1 DO DATBUF[I]←0;
04000 END;
04100 BPT←POINT(BPS,DATBUF[0],-1);
04200 FOR K←1 STEP SAMINT UNTIL 3*DATSIZ%(2*N) DO
04300 BEGIN
04400 SEGC←SEGC+1;
04500 DPYSET(DPYBUF);
04600 IF CHAR1="N" THEN OUTSTR('15&'12&CVS(SEGC))
04700 ELSE
04800 BEGIN
04900 AIVECT(-511,-500);
05000 DPYSST(DATIME&" Data file "&FILEI&" Segment "&CVS(SEGC)&" of "&cvs(segtot)
05100 &" SPACING "&CVS(SAMINT)&" M="&CVS(M));
05200 END;
05300 TIMSET;
05400 UNPACK(A,B,DATBUF,BPT,N,SAMINT);
05500 TIMDPY(100,-80,RTIM," Unpack time = ");
05600 IF CHAR1="Y" THEN
05700 BEGIN
05750 SETFORMAT(2,1);
05800 DUBDIS(A,B,N-1,-511,256,"Original segment");
05850 AIVECT(0,420);
05875 DPYSST("Segment length = "&CVF(N*SAMINT/20));
05887 IF CHAR2≠"Y" THEN
05900 DPYOUT(1);
06000 END;
06100 IF CHAR3="Y" THEN DPYDMP;
06150 IF SEGC<SEGCNT THEN
06175 BEGIN
06187 WAIT(2);
06193 IF INCHRS=" " THEN SEGCNT←0;
06196 END;
06200 IF CHAR2="Y" THEN IF SEGC≥SEGCNT THEN
06300 BEGIN
06320 TIMSET;
06340 WINDOW(A,B,N,3*N%32);
06360 TIMDPY(100,-40,RTIM,"WINDOW TIME = ");
06380 IF CEPS≠"Y" THEN DUBDIS(A,B,N-1,-511,0,"WINDOWED SEGMENT");
06400 TIMSET;
06500 RFOUR(A,B,M,0);
06600 TIMDPY(100,-110,RTIM," FFT time = ");
06700 TIMSET;
06800 POWER(A,B,C,N);
06900 TIMDPY(100,-140,RTIM," P.Spectrum time = ");
07000 IF CHAR5="Y" THEN ARRYOUT(AOFORM,C[0],256);
07100 IF CHAR1="Y" THEN
07200 BEGIN
07210 IF LOGFRE="Y" THEN
07220 LOGDIS(C,D,N,-511,-400,"Power spectrum against log frequency")
07230 ELSE
07300 ARRDIS(C,N,-511,-400,"Power spectrum with each major division = "&cvs(1016%samint)&" cycles per second");
07320 IF CEPS="Y" THEN
07325 BEGIN
07330 LOGMAG(A,B,M);
07335 RFOUR(A,B,M,1);
07340 DUBDIS(A,B,N-1,-511,0,"CEPSTRUM");
07345 END;
07350 DPYOUT(1);
07400 IF CHAR4="Y" THEN
07500 BEGIN
07600 AIVECT(0,-180);
07700 DPYSST("Space bar to continue");
07750 DPYOUT(1);
07800 IF(INCHRW)="S" THEN SEGCNT←CVD(INCHWL);
07900 END;
08100 END;
08200 IF CHAR3="Y" THEN DPYDMP;
08300 END;
08400 END;
08500 END;
08600
08700 DPYCLR;
08800 CLOSIN(AIFORM); CLOSO(AOFORM);
08900 RELEASE(AIFORM); RELEASE(AOFORM);
09000 END "FAST";